﻿using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using WpfAlgorithm1;

namespace Aggregation_Clustering
{
    class ClusterList
    {
        public ClusterList()
        {
            clusterList = new ArrayList();
            elementList = new ArrayList();
            //set the dummyCluster
            clusterList.Add(new Cluster(0));
        }
        //初始化所有S个Cluster以及初始值
        public void InitClusters(int s, int ubound, int lbound, double a1 = 1, double a2 = 0)
        {
            numS = s;

            UBound = ubound;
            LBound = new int[numS];

            A1 = a1;
            A2 = a2;
            VarianceConcerned = false;

            for (int i = 0; i < numS; i++)
            {
                LBound[i] = lbound;
            }
            Cluster DCluster = (Cluster)clusterList[0];
            numT = DCluster.num;//number of the elements in dummyCluster

            MoveVal = new double[numT + 1, numS + 1];//注意其实有s个cluster' 和1个Dcluster
            Move_Vari_Val = new double[numT + 1, numS + 1];//注意其实有s个cluster' 和1个Dcluster

            CandidateList = new ArrayList[numS + 1];
            for (int i = 0; i < numS + 1; i++)
            {
                CandidateList[i] = new ArrayList();
            }
            AddVal = new double[numT + 1, numS + 1];//注意其实有s个cluster' 和1个Dcluster
            Sum_Value = new double[numS + 1];
            Vari_Value = new double[numS + 1];


            TabuEnd = new double[numT + 1, numS + 1];//TS中的数组
            DropPenalty = new double[numS + 1];
            AddPenalty = new double[numS + 1];


            for (int i = 1; i <= numS; i++)
            {
                clusterList.Add(new Cluster(i));
            }
            ObjVal = 0;//将最初的cluster价值设为0
            //StartWeight = 0.7;//初始权重(为凸显某个因素的权)设为0.7
            // Weight = StartWeight;//weight不用求
            //printf “Initial Threshold Test Weight = ” StartWeight // Foregoing print is to report the initial value of Weight used in the Threshold Test 
            //large = double.MaxValue;

            //set elements' ID in initialization


        }

        //add Element to the dummy Cluster
        public void AddElement2defaultCluster(Element t)
        {
            Cluster tempcluster = (Cluster)clusterList[0];
            tempcluster.Add(t);
            elementList.Add(t);
            t.ElementID = elementList.IndexOf(t) + 1;
        }

        //4.3 First create 2-element clusters.
        public void Create2EleClusters()
        {
            //calculate the CheckBound
            Cluster DCluster = (Cluster)clusterList[0];
            numT = DCluster.num;//number of the elements in dummyCluster
            int CheckBound = numT / numS;
            int newLBound = (numT / numS) / 2;
            if (LBound[0] < newLBound)
            {
                for (int i = 0; i < numS; i++)
                {
                    LBound[i] = newLBound;//////////////////////////////////////////////////////////////////////////////////
                }
            }
            if (CheckBound >= 7)
                CheckBound = CheckBound - 2;
            CheckBound = MathHelper.Min(CheckBound, 10);//每个cluster中最多有10个元素
            ///////////////////////////////////////////////////////////////////
            ////////////////////////////////////////////////////////////////////////////////

            Element elementT0Best = null, elementT1Best = null, elementT2Best = null, elementT3Best = null;
            Element elementT0 = null, elementT1 = null, elementT2 = null, elementT3 = null;
            int TLeft;

            #region 填装所有CLuster
            for (int s = 1; s <= numS; s++)
            {
                double BestSumScore = double.MaxValue;

                elementT0 = DCluster.First;
                #region//Endwhile (t0 > 0)//当此次循环结束时当前Cluster（s）所有点的最优组合（t1,t2,t3,t0）就已找到。
                while (elementT0 != null)
                {
                    double Best1 = double.MaxValue;
                    double Best2 = double.MaxValue;
                    double Best3 = double.MaxValue;
                    elementT1 = null;
                    elementT2 = null;
                    elementT3 = null;
                    // for each t0, find the three best scores Score(t0,t1), Score(t0,t2) and Score(t0,t3). 

                    Element elementTCheck = DCluster.First;


                    #region//当这个while循环结束后对于固定t0已经找到了最小的
                    while (elementTCheck != null)
                    {
                        if (elementTCheck != elementT0)
                        {
                            double TrialVal = CalculateScore(elementT0, elementTCheck);
                            if (TrialVal < Best3)
                            {
                                if (TrialVal < Best2)
                                {
                                    Best3 = Best2;
                                    elementT3 = elementT2;
                                    if (TrialVal < Best1)
                                    {
                                        Best2 = Best1;
                                        elementT2 = elementT1;
                                        Best1 = TrialVal;
                                        elementT1 = elementTCheck;
                                    }
                                    else
                                    {
                                        Best2 = TrialVal;
                                        elementT2 = elementTCheck;
                                    }
                                }
                                else
                                {
                                    Best3 = TrialVal;
                                    elementT3 = elementTCheck;
                                }
                            }
                        }
                        elementTCheck = elementTCheck.After;
                    }
                    #endregion

                    if (elementT3 == null)
                    {
                        Best3 = 0;
                    }
                    if (elementT2 == null)
                    {
                        Best2 = 0;
                    }
                    if (elementT1 == null)
                    {
                        Best1 = 0;
                    }

                    double SumScore = Best1 + Best2 + Best3;// Score(t0,t1), Score(t0,t2) ,Score(t0,t3)
                    if (SumScore < BestSumScore)
                    {//若此时的SumScore小于之前的最优和
                        BestSumScore = SumScore; //则用新的结果更新最优和
                        elementT0Best = elementT0;
                        elementT1Best = elementT1;
                        elementT2Best = elementT2;
                        elementT3Best = elementT3;//记录t0的最优组合
                    }

                    elementT0 = elementT0.After;
                }
                #endregion
                double Try1 = double.MaxValue;
                double Try2 = double.MaxValue;
                double Try3 = double.MaxValue;
                //only has t0 left
                if (elementT1Best != null && elementT2Best != null)
                {
                    Try1 = CalculateScore(elementT0Best, elementT1Best) + CalculateScore(elementT0Best, elementT2Best) + CalculateScore(elementT1Best, elementT2Best);
                }
                if (elementT1Best != null && elementT3Best != null)
                {
                    Try2 = CalculateScore(elementT0Best, elementT1Best) + CalculateScore(elementT0Best, elementT3Best) + CalculateScore(elementT1Best, elementT3Best);
                }
                if (elementT2Best != null && elementT3Best != null)
                {
                    Try3 = CalculateScore(elementT0Best, elementT2Best) + CalculateScore(elementT0Best, elementT3Best) + CalculateScore(elementT2Best, elementT3Best);
                }
                #region 找出当前t0下最优的t1和t2
                elementT0 = elementT0Best;
                if (Try1 <= Try2 && Try1 <= Try3)//try1 is the best
                {
                    elementT1 = elementT1Best;
                    elementT2 = elementT2Best;
                }
                else
                {
                    if (Try2 <= Try1 && Try2 <= Try3)//try2 is the best
                    {
                        elementT1 = elementT1Best;
                        elementT2 = elementT3Best;
                    }
                    else
                    {
                        if (Try3 <= Try1 && Try3 <= Try2)//try3 is the best
                        {
                            elementT1 = elementT2Best;
                            elementT2 = elementT3Best;
                        }
                    }
                }
                // Now the pair t0,t1 creates the next 2-element Cluster(s)
                // Remove t0 and t1 from the dummy Cluster(0)
                #endregion

                #region 找出T1、T2中相对更优的T1
                if (CalculateScore(elementT0, elementT1) > CalculateScore(elementT0, elementT2))
                {
                    elementT1 = elementT2;
                }
                #endregion

                //
                // Now the pair t0,t1 creates the next 2-element Cluster(s)
                // Remove t0 and t1 from the dummy Cluster(0)
                elementT0.MoveTo(this, s);
                elementT1.MoveTo(this, s);

                //define the 2 unused parameters 
                //double SumMoveVal = 0.0;
                //double BestMoveVal = double.MaxValue;

                Element t = DCluster.First;
                while (t != null)
                {
                    AddVal[t.ElementID, s] = CalculateScore(t, elementT0) + CalculateScore(t, elementT1);
                    t = t.After;
                }

                #region Fill rest of Cluster s up to size CheckBound
                if (s <= numS)
                {
                    Cluster sCluster = getClusterbyID(s);
                    #region 该循环将当前Cluster填装至CheckBound
                    while (sCluster.num < CheckBound)
                    {
                        double BestVal = double.MaxValue;
                        Element tBest = null;
                        Element tCheck = DCluster.First;
                        double TrialVal;
                        while (tCheck != null)
                        {
                            TrialVal = AddVal[tCheck.ElementID, s];
                            if (TrialVal < BestVal)
                            {
                                BestVal = TrialVal;
                                tBest = tCheck;
                            }
                            tCheck = tCheck.After;
                        }
                        // Now add tBest to Cluster(s) by first dropping from Cluster(s0) 
                        tBest.MoveTo(this, s);

                        // Now update AddVal array.
                        t = DCluster.First;
                        while (t != null)
                        {
                            AddVal[t.ElementID, s] = AddVal[t.ElementID, s] + CalculateScore(t, tBest);
                            t = t.After;
                        }
                    }
                    #endregion
                }
                #endregion
            }
            #endregion
            //所有Cluster填装完毕
            //如果CheckBound大于2 则将其他元素回归

            #region 回归大于2的元素
            if (CheckBound > 2)
            {
                Element t0First, tsFirst, t0Last;//这些值在算法中用于整段的移动Element链到DCluster中
                //程序中用不到
                for (int s = 1; s <= numS; s++)
                {
                    Cluster sCluster = getClusterbyID(s);
                    t0First = sCluster.First;
                    for (int CheckSet = 0; CheckSet < CheckBound - 2; CheckSet++)
                    {
                        Element t = sCluster.First;
                        t.MoveTo(this, 0);
                    }
                    tsFirst = sCluster.First;
                    t0Last = tsFirst.Before;

                }
            }
            #endregion

            //     ObjVal = 0.0 from Initialization
            // 	Create arrays for each 2-element Cluster s
            for (int s = 1; s <= numS; s++)
            {
                Cluster sCluster = getClusterbyID(s);
                Element t0 = sCluster.First;
                Element t1 = t0.After;
                double sValue = CalculateScore(t0, t1);
                Sum_Value[s] = sValue;

                ObjVal += sValue;////////////ObjVal is defined in Chapter2   
                ///////////////////////////////////////////////////////////ObjVal = ∑ (FullValue(s): s = 1 to NumS)

                sCluster.num = 2;
                sCluster.NumberLinks = 1;

                //遍历求得两个矩阵值
                for (int i = 0; i <= numS; i++)
                {
                    Cluster tempCluster = getClusterbyID(i);
                    Element tempElement = tempCluster.First;
                    while (tempElement != null)
                    {
                        AddVal[tempElement.ElementID, s] = CalculateScore(tempElement, t0) + CalculateScore(tempElement, t1);
                        MoveVal[tempElement.ElementID, s] = (sValue + AddVal[tempElement.ElementID, s]) / 3.0 - sValue;//在下面会用到
                        tempElement = tempElement.After;
                    }

                }
            }
            TLeft = numT - 2 * numS;
            try
            {
                CompleteTheInit(TLeft);
            }
            catch (Exception e)
            {
                throw e;
            }
        }

        public void CompleteTheInit(int tleft)
        {

            ObjValList = new ArrayList();
            Obj_Vari_ValList = new ArrayList();


            ObjValList.Clear();
            Obj_Vari_ValList.Clear();


            int CheckBound;
            if (tleft != numT - 2 * numS)
            {
                throw new Exception("ERROR: TLEFT is wrong." +
                    " TLEFT =" + tleft.ToString() +
                    " numT =" + numT.ToString() +
                    " numS =" + numS.ToString());
                //                Print “Error: TLeft is wrong”
                //Print “TLeft = ” TLeft
                //Print “NumT = ” NumT
                //Print “NumS = ” NumS”
            }//首先检查TLeft的值是否正确，错误则报错并退出。

            int BndLeft = 0;
            for (int i = 0; i < numS; i++)
            {
                BndLeft += (LBound[i] - 2);
            }
            int NumLeft = tleft;

            Cluster DCluster = getClusterbyID(0);

            #region// Each iteration assigns an element t0 of Cluster(0) to some Cluster(s) by Min Regret.
            for (int i = 0; i < tleft; i++)
            {
                if (BndLeft == NumLeft)
                {
                    CheckBound = LBound[0];//warning:
                    //each Cluster has the same LBound at current Time
                    /////////////////////////////////////////////////////////////////////////////////////////////////
                }
                else
                {
                    CheckBound = UBound;
                }
                double MinRegret = double.MaxValue;
                Element tCheck = DCluster.First;
                Element tBest = tCheck; ;
                int sBest = 0;

                #region//找到做小regret MinRegret 的element tBest和对应的cluster sBest
                while (tCheck != null)
                {
                    double BestVal1 = double.MaxValue;
                    double BestVal2 = double.MaxValue;
                    int s1Best = 0;
                    int s2Best = 0;
                    for (int s = 1; s <= numS; s++)
                    {
                        Cluster curCluster = getClusterbyID(s);
                        if (curCluster.num < CheckBound)
                        {
                            double TrialVal = MoveVal[tCheck.ElementID, s];
                            if (TrialVal < BestVal2)
                            {
                                if (TrialVal < BestVal1)
                                {
                                    BestVal2 = BestVal1;
                                    s2Best = s1Best;
                                    BestVal1 = TrialVal;
                                    s1Best = s;//replace s2 with s1 ,s1 with s
                                }
                                else
                                {
                                    BestVal2 = TrialVal;
                                    s2Best = s;//replace s2with s
                                }
                            }
                        }
                    }//find the best 2 Cluster s1best and s2best

                    double Regret = BestVal1 - BestVal2;
                    if (Regret < MinRegret)
                    {
                        MinRegret = Regret;
                        tBest = tCheck;
                        sBest = s1Best;
                    }
                    tCheck = tCheck.After;
                }
                #endregion

                //MOve tBest from cLuster(0) to Cluster(sBest)
                tBest.MoveTo(this, sBest);
                Sum_Value[sBest] += AddVal[tBest.ElementID, sBest];
                double DeltaVal = MoveVal[tBest.ElementID, sBest];

                ObjVal += DeltaVal;////////////ObjVal is defined in Chapter2   and assigned in Chapter4.3
                ///////////////////////////////////////////////////////////ObjVal = ∑ (FullValue(s): s = 1 to NumS)
                /////////////////////////////////////////////////ObjVal = ObjVal + sValue

                Cluster Sa = getClusterbyID(sBest);
                double SAValue = Sum_Value[sBest];
                double PresentFullValueA = SAValue / Sa.NumberLinks;
                //update values
                for (int t = 1; t <= numT; t++)
                {
                    Element elementT = getElementbyID(t);
                    AddVal[t, sBest] += CalculateScore(elementT, tBest);
                    double FullAddValA = ((SAValue + AddVal[t, sBest]) / Sa.NumberLinks) - PresentFullValueA;
                    MoveVal[t, sBest] = FullAddValA;
                }//update done
                // Moreover, we only move elements from the dummy cluster, so remaining 
                // MoveVal(t,sA) values are irrelevant until all elements are finally assigned.

                NumLeft = NumLeft - 1;
                if (Sa.num <= LBound[0])
                {
                    BndLeft--;
                }
            }
            #endregion

            if (DCluster.First != null)
            {
                throw new Exception("Error: First(0) ≠ 0:  First(0) = " + DCluster.First.ElementID.ToString());
            }

        }

        private void InitCandidateList()
        //initial CandidateList by All Random Strategy
        //pick 20% availables moves
        {
            Random ra = new Random();
            int availabelMoveNum=numS*numT/5;
            for (int i = 0; i < availabelMoveNum; i++)
            {
                int toclusterID= ra.Next(1, numS + 1);
                Candidate tempcandidate = new Candidate(ra.Next(1, numT + 1),toclusterID);
                CandidateList[toclusterID].Add(tempcandidate);                  
            }
        }
        public void InitTSPhase()
        {
            for (int t = 1; t <= numT; t++)
            {
                for (int s = 1; s <= numS; s++)
                {
                    AddVal[t, s] = 0;
                    TabuEnd[t, s] = 0;
                }
            }//reset addval& tabuend
            for (int t = 1; t <= numT; t++)
            {
                int sA = getElementbyID(t).ClusterID;
                for (int t0 = 1; t0 <= numT; t0++)
                {
                    AddVal[t0, sA] += CalculateScore(getElementbyID(t0), getElementbyID(t));
                }
            }//determine the AddVal

            ObjVal = 0;
            Obj_Vari_Val = 0;
            for (int s = 1; s <= numS; s++)
            {
                Cluster curCluster = getClusterbyID(s);
                int CheckLinkS = 0;
                double ValS = 0;
                Element t1 = curCluster.First;
                Element t2 = null;
                while (t1 != null)
                {
                    t2 = t1.After;
                    while (t2 != null)
                    {
                        CheckLinkS++;
                        ValS += CalculateScore(t1, t2);
                        t2 = t2.After;
                    }
                    t1 = t1.After;
                }
                Sum_Value[s] = ValS;
                curCluster.NumberLinks = CheckLinkS;
                /// Here comes BugCheck1
                /////////////////////////////////////////////////////////////////////////////
                int curNumLinks = curCluster.num * (curCluster.num - 1) / 2;
                if (CheckLinkS != curNumLinks)
                {
                    Exception e = new Exception("BugCheck1: Cluster NO." + s.ToString() + " CheckLinkS:" + CheckLinkS + "; ElementNum:" + curCluster.num);
                    throw e;
                }
                double FullValS = ValS / CheckLinkS;

                Vari_Value[s] = VariValuesCalculator(curCluster);
                Obj_Vari_Val += Vari_Value[s];

                if (!VarianceConcerned)
                {
                    ObjVal += FullValS;
                }
                else
                {
                    ObjVal += FullValS * A1 + Vari_Value[s] * A2;
                }
            }
            //initial CandidateList by All Random Strategy
            //pick 20% availables moves
            InitCandidateList();


            DropPenVal = 0.5 * double.MaxValue;
            #region //here Prepare info(MoveValue Move_Vari_Val) for s1 to evaluate adding some element t2 to it.
            //value pairs (t2,s1) are in CandidateList[s1]
            for (int s1 = 1; s1 <= numS; s1++)
            {
                Cluster clusterS1 = getClusterbyID(s1);
                int Link1 = clusterS1.NumberLinks;
                int NewLink1 = Link1 + clusterS1.num;
                double Val1 = Sum_Value[s1];
                double FullVal1 = Val1 / Link1;
                foreach (Candidate candidate in CandidateList[s1])
                {
                    int t2 = candidate.EleID;
                    Element elementT2 = getElementbyID(t2);
                    int s2 = elementT2.ClusterID;
                    Cluster clusterS2 = getClusterbyID(s2);
                    if (s1 == s2)
                    {
                        MoveVal[elementT2.ElementID, s1] = 0;
                        Move_Vari_Val[elementT2.ElementID, s1] = 0;
                    }
                    else if (clusterS2.num == 2)//set penalty if only has 2 elements
                    {
                        MoveVal[elementT2.ElementID, s1] = DropPenVal;
                        Move_Vari_Val[elementT2.ElementID, s1] = DropPenVal;
                    }
                    else// set new Value
                    {
                        int Link2 = clusterS2.NumberLinks;
                        int NewLink2 = Link2 - (clusterS2.num - 1);
                        double Val2 = Sum_Value[s2];
                        double FullVal2 = Val2 / Link2;
                        
                        double NewFullVal1 = (Val1 + AddVal[t2, s1]) / NewLink1;
                        double DeltaAdd1 = NewFullVal1 - FullVal1;
                        
                        double NewFullVal2 = (Val2 - AddVal[t2, s2]) / NewLink2;
                        double DeltaDrop2 = NewFullVal2 - FullVal2;
                        
                        MoveVal[t2, s1] = DeltaAdd1 + DeltaDrop2;
                        Move_Vari_Val[t2, s1] = Move_Vari_ValueCalculator(elementT2, clusterS1);
                    }
                }      
            }

            #region //previous version
            //for (int s1 = 1; s1 <= numS; s1++)
            //{
            //    Cluster clusterS1 = getClusterbyID(s1);
            //    int Link1 = clusterS1.NumberLinks;
            //    int NewLink1 = Link1 + clusterS1.num;
            //    double Val1 = Sum_Value[s1];
            //    double FullVal1 = Val1 / Link1;

            //    for (int s2 = 1; s2 <= numS; s2++)
            //    {
            //        Cluster clusterS2 = getClusterbyID(s2);
            //        if (s2 == s1)
            //        {
            //            Element elementT2 = clusterS2.First;
            //            while (elementT2 != null)
            //            {
            //                MoveVal[elementT2.ElementID, s1] = 0;
            //                Move_Vari_Val[elementT2.ElementID, s1] = 0;
            //                elementT2 = elementT2.After;
            //            }
            //        }
            //        else
            //        {
            //            if (clusterS2.num == 2)//set penalty if only has 2 elements
            //            {
            //                Element elementT2 = clusterS2.First;
            //                while (elementT2 != null)
            //                {
            //                    MoveVal[elementT2.ElementID, s1] = DropPenVal;
            //                    Move_Vari_Val[elementT2.ElementID, s1] = DropPenVal;
            //                    elementT2 = elementT2.After;
            //                }
            //            }
            //            else//or set the new moveval
            //            {
            //                int Link2 = clusterS2.NumberLinks;
            //                int NewLink2 = Link2 - (clusterS2.num - 1);
            //                double Val2 = Sum_Value[s2];
            //                double FullVal2 = Val2 / Link2;

            //                Element elementT2 = clusterS2.First;
            //                while (elementT2 != null)
            //                {
            //                    double NewFullVal1 = (Val1 + AddVal[elementT2.ElementID, s1]) / NewLink1;
            //                    double DeltaAdd1 = NewFullVal1 - FullVal1;

            //                    double NewFullVal2 = (Val2 - AddVal[elementT2.ElementID, s2]) / NewLink2;
            //                    double DeltaDrop2 = NewFullVal2 - FullVal2;

            //                    MoveVal[elementT2.ElementID, s1] = DeltaAdd1 + DeltaDrop2;
            //                    Move_Vari_Val[elementT2.ElementID, s1] = Move_Vari_ValueCalculator(elementT2, clusterS1);
            //                    elementT2 = elementT2.After;
            //                }
            //            }
            //        }
            //    }
            //}
            #endregion
            #endregion

            #region 4.7# Prepare Penalties and Inducements for Tabu Search
            AddPenVal = 0.5 * double.MaxValue;
            DropPenVal = 0.5 * double.MaxValue;

            AddInduceVal = -double.MaxValue / 5;
            DropInduceVal = -double.MaxValue / 5;

            NumUBViolate = 0;
            NumLBViolate = 0;

            for (int s = 1; s <= numS; s++)
            {
                Cluster clusterS = getClusterbyID(s);
                if (clusterS.num >= UBound)
                {
                    int AddViolation1 = 1 + clusterS.num - UBound;
                    ///////////////////////////////////
                    //////////////////NumUBViolate??????????????????????????????????????????????????????????????????????????
                    //NumUBViolate += AddViolation1;
                    AddPenalty[s] = AddPenVal * AddViolation1;
                    if (clusterS.num > UBound)
                    {
                        int AddViolation0 = clusterS.num - UBound;
                        NumUBViolate += AddViolation0;
                        DropPenalty[s] = DropInduceVal * AddViolation0;
                    }
                    else
                    {
                        DropPenalty[s] = 0;
                    }
                }
                else if (clusterS.num <= LBound[s - 1])
                {
                    int DropViolation1 = 1 + LBound[s - 1] - clusterS.num;
                    DropPenalty[s] = DropPenVal * DropViolation1;
                    if (clusterS.num < LBound[s - 1])
                    {
                        int DropViolation0 = LBound[s - 1] - clusterS.num;
                        NumLBViolate += DropViolation0;
                        AddPenalty[s] = AddInduceVal * DropViolation0;
                    }
                    else
                    {
                        AddPenalty[s] = 0;
                    }
                }
                else
                {
                    AddPenalty[s] = 0;
                    DropPenalty[s] = 0;
                }
            }
            NumAllViolate = NumLBViolate + NumUBViolate;
            #endregion

            #region 4.8# prepare parameters for TS

            SmallTabuTenure = MathHelper.Min(3, numS / 3);
            SmallTabuTenure = MathHelper.Max(SmallTabuTenure, 1);

            double TabuTenureLim = 1 + 0.4 * (numT - SmallTabuTenure) * (numS - 1);
            double TabuTenurePref = 1 + 0.25 * (numT - SmallTabuTenure) * (numS - 1);

            TabuTenure = MathHelper.Max(TabuTenurePref, 15);
            TabuTenure = MathHelper.Min(TabuTenure, TabuTenureLim);

            #endregion
        }


        // #4.80 Launch Main Body of Tabu Search Improvement Method
        public void ExecTSPhase(int MaxIter = 1000, int NoImproveLimit = 100)
        // Set limits on the number of iterations and of non-improving iterations by default value 1000,100
        // For extreme:
        // MaxIter = 50000
        // NoImproveLimit = 200//////
        {
            MoveLogList = new ArrayList();
            ClusterLogList = new ArrayList();

            BestImprove = true;
            LastImproveIter = 0;
            Iter = 1;
            BestObjVal = ObjVal;
            BestVariVal = Obj_Vari_Val;
            NumNewBest = 0;
            LastBestIter = 0;
            FirstBestIter = -1;

            double BestMoveVal, Aspiration;
            Cluster clusterS0, clusterSA, clusterSBest = null;
            Element elementT0, elementTBest = null;
            int[] bestClusterID = new int[numT];
            int TabuIter;
            while (Iter < MaxIter)//main loop for ts  
            {
                if ((Iter - LastImproveIter) < NoImproveLimit)
                {
                    #region//////improvement phases
                    if (BestImprove)
                    {
                        TabuIter = int.MaxValue;
                    }
                    else
                    {
                        TabuIter = Iter;
                    }
                    // Identify the best non-tabu move 
                    BestMoveVal = double.MaxValue;
                    Aspiration = BestObjVal - ObjVal;

                    #region #4.80 find the best elementT and ClusterS
                    bool Drop0Feas, Drop1Feas;
                    for (int s = 1; s <= numS; s++)
                    {
                        Cluster clusterS = getClusterbyID(s);
                        foreach (Candidate candidate in CandidateList[s])
                        {
                            int T0=candidate.EleID;
                            elementT0=getElementbyID(T0);
                            int S0=elementT0.ClusterID;
                            clusterS0 = getClusterbyID(S0);

                            if (clusterS0.num>2)
                            {
                                double s0DropPen = DropPenalty[S0];
                                if (TabuEnd[T0, S0] < TabuIter)
                                {
                                    if (s != S0)
                                    {
                                        double TrueMoveVal = MoveVal[elementT0.ElementID, s];
                                        double TrueMove_Vari_Val = Move_Vari_Val[elementT0.ElementID, s];
                                        double TestMoveVal;
                                        if (VarianceConcerned)
                                        {
                                            TrueMoveVal = A2 * TrueMove_Vari_Val + A1 * TrueMoveVal;
                                        }
                                        TestMoveVal = TrueMoveVal + AddPenalty[s] + s0DropPen;

                                        if (TestMoveVal < BestMoveVal)
                                        {
                                            Drop0Feas = false;
                                            Drop1Feas = false;
                                            //if element can be moved from S0
                                            if (clusterS0.num > LBound[S0 - 1])
                                            {
                                                if (NumAllViolate == 0
                                                    ||
                                                    (NumAllViolate == 1 && clusterS0.num == UBound + 1))
                                                {
                                                    Drop0Feas = true;
                                                }
                                                else
                                                {
                                                    if (NumAllViolate <= 2)
                                                    {
                                                        if (NumLBViolate == 1)
                                                        {
                                                            if (NumAllViolate == 1 || clusterS0.num == UBound + 1)
                                                            {
                                                                Drop1Feas = true;
                                                            }
                                                        }
                                                    }
                                                }
                                            }
                                            if (TabuEnd[T0, s] < TabuIter
                                                || ((TrueMoveVal < Aspiration)
                                                         && ((Drop0Feas == true && clusterS.num < UBound)
                                                                ||
                                                                (Drop1Feas == true && clusterS.num == LBound[s - 1] - 1)
                                                            )
                                                    )
                                                )
                                            {
                                                BestMoveVal = TestMoveVal;
                                                elementTBest = elementT0;
                                                clusterSBest = clusterS;
                                            }
                                        }
                                    }
                                } 
                            }
                        }
                    }
                    #region previous Version
                    //for (int s0 = 1; s0 <= numS; s0++)
                            //{
                            //    clusterS0 = getClusterbyID(s0);
                            //    if (clusterS0.num > 2)
                            //    {
                            //        Drop0Feas = false;
                            //        Drop1Feas = false;
                            //        //if element can be moved from S0
                            //        if (clusterS0.num > LBound[s0 - 1])
                            //        {
                            //            if (NumAllViolate == 0
                            //                ||
                            //                (NumAllViolate == 1 && clusterS0.num == UBound + 1))
                            //            {
                            //                Drop0Feas = true;
                            //            }
                            //            else
                            //            {
                            //                if (NumAllViolate <= 2)
                            //                {
                            //                    if (NumLBViolate == 1)
                            //                    {
                            //                        if (NumAllViolate == 1 || clusterS0.num == UBound + 1)
                            //                        {
                            //                            Drop1Feas = true;
                            //                        }
                            //                    }
                            //                }
                            //            }
                            //        }
                            //        double s0DropPen = DropPenalty[s0];

                            //        elementT0 = clusterS0.First;
                            //        while (elementT0 != null)
                            //        {
                            //            if (TabuEnd[elementT0.ElementID, s0] < TabuIter)
                            //            {
                            //                for (int s = 1; s <= numS; s++)
                            //                {
                            //                    Cluster clusterS = getClusterbyID(s);
                            //                    if (s != s0)
                            //                    {
                            //                        double TrueMoveVal = MoveVal[elementT0.ElementID, s];
                            //                        double TrueMove_Vari_Val = Move_Vari_Val[elementT0.ElementID, s];
                            //                        double TestMoveVal;
                            //                        if (VarianceConcerned)
                            //                        {
                            //                            TrueMoveVal = A2 * TrueMove_Vari_Val + A1 * TrueMoveVal;
                            //                        }
                            //                        TestMoveVal = TrueMoveVal + AddPenalty[s] + s0DropPen;

                            //                        if (TestMoveVal < BestMoveVal)
                            //                        {
                            //                            if (TabuEnd[elementT0.ElementID, s] < TabuIter
                            //                                || ((TrueMoveVal < Aspiration)
                            //                                   && ((Drop0Feas == true && clusterS.num < UBound)
                            //                                        ||
                            //                                       (Drop1Feas == true && clusterS.num == LBound[s - 1] - 1)
                            //                                     )
                            //                                  )
                            //                              )
                            //                            {
                            //                                BestMoveVal = TestMoveVal;
                            //                                elementTBest = elementT0;
                            //                                clusterSBest = clusterS;
                            //                            }
                            //                        }
                            //                    }
                            //                }
                            //            }
                            //            elementT0 = elementT0.After;
                            //        }
                            //    }
                    //}
                    #endregion
                    #endregion

                    #region #4.9 Move and Update Sum_Values
                    //source
                    elementT0 = elementTBest;
                    clusterS0 = getClusterbyID(elementT0.ClusterID);
                    //target
                    clusterSA = clusterSBest;

                    elementT0.MoveTo(this, clusterSA.ID);
                    Sum_Value[clusterSA.ID] += AddVal[elementT0.ElementID, clusterSA.ID];
                    Sum_Value[clusterS0.ID] -= AddVal[elementT0.ElementID, clusterS0.ID];

                    Vari_Value[clusterSA.ID] = VariValuesCalculator(clusterSA);
                    Vari_Value[clusterS0.ID] = VariValuesCalculator(clusterS0);

                    #region//here comes Bug Check2
                    if (clusterSA.NumberLinks != clusterSA.num * (clusterSA.num - 1) / 2 || clusterS0.NumberLinks != clusterS0.num * (clusterS0.num - 1) / 2)
                    {
                        Exception e = new Exception("BugCheck2:"
                        + "\r\nCluster NO." + clusterSA.ID.ToString() + " CheckLinkS:" + clusterSA.NumberLinks + "; ElementNum:" + clusterSA.num
                        + "\r\nCluster NO." + clusterS0.ID.ToString() + " CheckLinkS:" + clusterS0.NumberLinks + "; ElementNum:" + clusterS0.num);
                        throw e;
                    }
                    #endregion

                    int Link0 = clusterS0.NumberLinks;
                    int LinkA = clusterSA.NumberLinks;
                    double ValA = Sum_Value[clusterSA.ID];
                    double Val0 = Sum_Value[clusterS0.ID];

                    double FullVal0 = Val0 / Link0;
                    int NewAddLink0 = Link0 + clusterS0.num;
                    int NewDropLink0 = Link0 - clusterS0.num + 1;

                    double FullValA = ValA / LinkA;
                    int NewAddLinkA = LinkA + clusterSA.num;
                    int NewDropLinkA = LinkA - clusterSA.num + 1;

                    #region//Update Violates
                    if (clusterS0.num >= UBound)
                    {
                        NumUBViolate--;
                        NumAllViolate--;
                    }
                    else if (clusterS0.num < LBound[clusterS0.ID - 1])
                    {
                        NumLBViolate++;
                        NumAllViolate++;
                    }
                    if (clusterSA.num > UBound)
                    {
                        NumUBViolate++;
                        NumAllViolate++;
                    }
                    else if (clusterSA.num <= LBound[clusterSA.ID - 1])
                    {
                        NumLBViolate--;
                        NumAllViolate--;
                    }
                    #endregion

                    double DeltaVal = MoveVal[elementT0.ElementID, clusterSA.ID];
                    double Delta_Vari_Val = Move_Vari_Val[elementT0.ElementID, clusterSA.ID];

                    if (VarianceConcerned)
                    {
                        DeltaVal = A1 * DeltaVal + A2 * Delta_Vari_Val;
                    }

                    ObjVal += DeltaVal;

                    Obj_Vari_Val += Delta_Vari_Val;

                    //// Store new Best solution if solution improves
                    if (NumAllViolate == 0)
                    {
                        if (DeltaVal < 0)
                        {
                            if (ObjVal <= BestObjVal)
                            {
                                LastImproveIter = Iter;
                                BestImprove = true;
                                BestObjVal = ObjVal;
                                BestVariVal = Obj_Vari_Val;
                            }
                            else
                            {
                                BestImprove = false;
                            }
                        }
                        else if (BestImprove == true)
                        {
                            //store the best solution
                            Element tempElement;
                            for (int t = 0; t < numT; t++)
                            {
                                tempElement = getElementbyID(t + 1);
                                bestClusterID[t] = tempElement.ClusterID;
                            }
                            bestClusterID[elementT0.ElementID - 1] = clusterS0.ID;

                            BestImprove = false;
                            NumNewBest++;
                            LastBestIter = Iter;
                            if (FirstBestIter < 0) FirstBestIter = Iter;
                        }
                    }
                    #endregion

                    #region #4.10 Compute New AddVal and MoveVal quantities for Next Tabu Iteration
                    //Update All AddVal and Candidates
                    for (int s2 = 1; s2 <= numS; s2++)
                    {
                        Cluster clusterS2 = getClusterbyID(s2);
                        Element elementT2 = clusterS2.First;
                        while (elementT2 != null)
                        {
                            AddVal[elementT2.ElementID, clusterS0.ID] -= CalculateScore(elementT2, elementT0);
                            AddVal[elementT2.ElementID, clusterSA.ID] += CalculateScore(elementT2, elementT0);
                            elementT2 = elementT2.After;
                        }
                    }
                    int[] clusternum = new int[numS];
                    for (int s = 1; s <= numS; s++)
                    {
                        Cluster clustertemp = getClusterbyID(s);
                        clusternum[s - 1] = clustertemp.num;
                        CandidateList[s].Clear();
                    }
                    for (int t2 = 1; t2 <= numT; t2++)
                    {
                        Element elementT2 = getElementbyID(t2);
                        int s2 = elementT2.ClusterID;
                        Cluster clusterS2 = getClusterbyID(s2);
                        int Link2 = clusterS2.NumberLinks;
                        int NewLink2 = Link2 - (clusterS2.num - 1);
                        double Val2 = Sum_Value[s2];
                        double FullVal2 = Val2 / Link2;

                        double[] tempfullVal = new double[numS];
                        double tempfullTotal = 0;
                        for (int s = 1; s <= numS; s++)
                        {
                            tempfullVal[s - 1] = AddVal[t2, s] / clusternum[s - 1];
                            if (s!=s2)
                            {
                                tempfullTotal += tempfullVal[s - 1]; 
                            }
                        }
                        tempfullTotal = tempfullTotal / (numS - 1);
                        for (int s = 1; s <= numS; s++)
                        {
                            if (s != s2)
                            {
                                if (tempfullVal[s - 1] < tempfullTotal)
                                {
                                    CandidateList[s].Add(new Candidate(t2,s));
                                    
                                    Cluster clusterS = getClusterbyID(s);
                                    int Link1 = clusterS.NumberLinks;
                                    int NewLink1 = Link1 + clusterS.num;
                                    double Val1 = Sum_Value[s];
                                    double FullVal1 = Val1 / Link1;


                                    double NewFullVal1 = (Val1 + AddVal[t2, s]) / NewLink1;
                                    double DeltaAdd1 = NewFullVal1 - FullVal1;

                                    double NewFullVal2 = (Val2 - AddVal[t2, s2]) / NewLink2;
                                    double DeltaDrop2 = NewFullVal2 - FullVal2;

                                    MoveVal[t2, s] = DeltaAdd1 + DeltaDrop2;
                                    Move_Vari_Val[t2, s] = Move_Vari_ValueCalculator(elementT2, clusterS);


                                }
                            }
                        }
                    }
                    #region previous Version
                    //for (int s2 = 1; s2 <= numS; s2++)
                    //{
                    //    Cluster clusterS2 = getClusterbyID(s2);
                    //    if (clusterS2.num == 2)
                    //    {
                    //        Element elementT2 = clusterS2.First;
                    //        while (elementT2 != null)
                    //        {
                    //            AddVal[elementT2.ElementID, clusterS0.ID] -= CalculateScore(elementT2, elementT0);
                    //            AddVal[elementT2.ElementID, clusterSA.ID] += CalculateScore(elementT2, elementT0);
                    //            MoveVal[elementT2.ElementID, clusterS0.ID] = DropPenVal;
                    //            MoveVal[elementT2.ElementID, clusterSA.ID] = DropPenVal;

                    //            Move_Vari_Val[elementT2.ElementID, clusterS0.ID] = DropPenVal;
                    //            Move_Vari_Val[elementT2.ElementID, clusterSA.ID] = DropPenVal;
                    //            elementT2 = elementT2.After;
                    //        }
                    //    }
                    //    else
                    //    {
                    //        int Link2 = clusterS2.NumberLinks;
                    //        int NewLink2 = Link2 - clusterS2.num + 1;
                    //        double Val2 = Sum_Value[s2];
                    //        double FullVal2 = Val2 / Link2;

                    //        Element elementT2 = clusterS2.First;
                    //        while (elementT2 != null)
                    //        {
                    //            AddVal[elementT2.ElementID, clusterS0.ID] -= CalculateScore(elementT2, elementT0);
                    //            AddVal[elementT2.ElementID, clusterSA.ID] += CalculateScore(elementT2, elementT0);

                    //            double NewFullVal2 = (Val2 - AddVal[elementT2.ElementID, s2]) / NewLink2;
                    //            double DeltaDrop2 = NewFullVal2 - FullVal2;

                    //            double NewFullVal0 = (Val0 + AddVal[elementT2.ElementID, clusterS0.ID]) / NewAddLink0;
                    //            double DeltaAdd0 = NewFullVal0 - FullVal0;
                    //            MoveVal[elementT2.ElementID, clusterS0.ID] = DeltaAdd0 + DeltaDrop2;
                    //            Move_Vari_Val[elementT2.ElementID, clusterS0.ID] = Move_Vari_ValueCalculator(elementT2, clusterS0);

                    //            double NewFullValA = (ValA + AddVal[elementT2.ElementID, clusterSA.ID]) / NewAddLinkA;
                    //            double DeltaAddA = NewFullValA - FullValA;
                    //            MoveVal[elementT2.ElementID, clusterSA.ID] = DeltaAddA + DeltaDrop2;
                    //            Move_Vari_Val[elementT2.ElementID, clusterSA.ID] = Move_Vari_ValueCalculator(elementT2, clusterSA);


                    //            elementT2 = elementT2.After;
                    //        }
                    //    }
                    //}
                    ////for elements from s0/sA
                    //for (int s1 = 1; s1 <= numS; s1++)
                    //{
                    //    Cluster clusterS1 = getClusterbyID(s1);
                    //    Element elementT2;
                    //    int Link1 = clusterS1.NumberLinks;
                    //    int NewLink1 = Link1 + clusterS1.num;
                    //    double Val1 = Sum_Value[s1];
                    //    double FullVal1 = Val1 / Link1;

                    //    #region// from s0 to s1

                    //    if (s1 == clusterS0.ID)
                    //    {
                    //        elementT2 = clusterS0.First;
                    //        while (elementT2 != null)
                    //        {
                    //            MoveVal[elementT2.ElementID, s1] = 0;
                    //            Move_Vari_Val[elementT2.ElementID, s1] = 0;
                    //            elementT2 = elementT2.After;
                    //        }
                    //    }
                    //    else
                    //    {
                    //        if (clusterS0.num == 2)
                    //        {
                    //            elementT2 = clusterS0.First;
                    //            while (elementT2 != null)
                    //            {
                    //                MoveVal[elementT2.ElementID, s1] = DropPenVal;
                    //                Move_Vari_Val[elementT2.ElementID, s1] = DropPenVal;

                    //                elementT2 = elementT2.After;
                    //            }
                    //        }
                    //        else
                    //        {
                    //            elementT2 = clusterS0.First;
                    //            while (elementT2 != null)
                    //            {
                    //                double NewFullVal0 = (Val0 - AddVal[elementT2.ElementID, clusterS0.ID]) / NewDropLink0;
                    //                double DeltaDrop0 = NewFullVal0 - FullVal0;

                    //                double NewFullVal1 = (Val1 + AddVal[elementT2.ElementID, clusterS1.ID]) / NewLink1;
                    //                double DeltaAdd1 = NewFullVal1 - FullVal1;
                    //                MoveVal[elementT2.ElementID, s1] = DeltaAdd1 + DeltaDrop0;
                    //                Move_Vari_Val[elementT2.ElementID, s1] = Move_Vari_ValueCalculator(elementT2, clusterS1);
                    //                elementT2 = elementT2.After;
                    //            }
                    //        }
                    //    }
                    //    #endregion

                    //    #region// from sA to s1
                    //    if (s1 == clusterSA.ID)
                    //    {
                    //        elementT2 = clusterSA.First;
                    //        while (elementT2 != null)
                    //        {
                    //            MoveVal[elementT2.ElementID, s1] = 0;
                    //            Move_Vari_Val[elementT2.ElementID, s1] = 0;
                    //            elementT2 = elementT2.After;
                    //        }
                    //    }
                    //    else
                    //    {
                    //        elementT2 = clusterSA.First;
                    //        while (elementT2 != null)
                    //        {
                    //            double NewFullValA = (ValA - AddVal[elementT2.ElementID, clusterSA.ID]) / NewDropLinkA;
                    //            double DeltadRropA = NewFullValA - FullValA;

                    //            double NewFullVal1 = (Val1 + AddVal[elementT2.ElementID, s1]) / NewLink1;
                    //            double DeltaAdd1 = NewFullVal1 - FullVal1;

                    //            MoveVal[elementT2.ElementID, s1] = DeltadRropA + DeltaAdd1;
                    //            Move_Vari_Val[elementT2.ElementID, s1] = Move_Vari_ValueCalculator(elementT2, clusterS1);

                    //            elementT2 = elementT2.After;
                    //        }
                    //    }
                    //    #endregion
                    //}

                    #endregion
                    string BugCheckLog;// = BugCheck3();
                    //if (BugCheckLog != "")
                    //{
                    //    BugCheckLog = "#Error in Iter No." + Iter.ToString() + "\r\n#While move Element(" + elementT0.ElementID.ToString() + ") from Cluster(" + clusterS0.ID.ToString() + ") to Cluster(" + clusterSA.ID.ToString() + ")\r\n" + BugCheckLog;
                    //}
                    //else
                    //{
                    BugCheckLog = "#Success in Iter No." + Iter.ToString() + "\r\n#While move Element(" + elementT0.ElementID.ToString() + ") from Cluster(" + clusterS0.ID.ToString() + ") to Cluster(" + clusterSA.ID.ToString() + ")" + "  MoveVal = " + DeltaVal.ToString("0.0000") + "  Move_Vari_Val = " + Delta_Vari_Val.ToString("0.0000") + "\r\n";
                    //}
                    BugCheckLog = BugCheckLog + "#Now ObjVal = " + ObjVal.ToString("0.0000");
                    if (VarianceConcerned)
                    {
                        BugCheckLog += "(VarianceConcerned)";
                    }
                    BugCheckLog += "\r\n";
                    ObjValList.Add(ObjVal);
                    Obj_Vari_ValList.Add(Obj_Vari_Val);
                    MoveLogList.Add(BugCheckLog);
                    #endregion

                    #region #4.11 Update Penalties and Inducements for sA and s0
                    if (clusterSA.num >= UBound)
                    {
                        int AddViolation1 = 1 + clusterSA.num - UBound;
                        AddPenalty[clusterSA.ID] = AddPenVal * AddViolation1;
                        if (clusterSA.num > UBound)
                        {
                            int AddViolation0 = clusterSA.num - UBound;
                            DropPenalty[clusterSA.ID] = DropInduceVal * AddViolation0;
                        }
                        else
                        {
                            DropPenalty[clusterSA.ID] = 0;
                        }
                    }
                    else if (clusterSA.num <= LBound[clusterSA.ID - 1])
                    {
                        int DropViolation1 = 1 + LBound[clusterSA.ID - 1] - clusterSA.num;
                        DropPenalty[clusterSA.ID] = DropPenVal * DropViolation1;
                        if (clusterSA.num < LBound[clusterSA.ID - 1])
                        {
                            int DropViolation0 = LBound[clusterSA.ID - 1] - clusterSA.num;
                            AddPenalty[clusterSA.ID] = AddInduceVal * DropViolation0;
                        }
                        else
                        {
                            AddPenalty[clusterSA.ID] = 0;
                        }
                    }
                    else
                    {
                        AddPenalty[clusterSA.ID] = 0;
                        DropPenalty[clusterSA.ID] = 0;
                    }

                    if (clusterS0.num >= UBound)
                    {
                        int AddViolation1 = 1 + clusterS0.num - UBound;
                        AddPenalty[clusterS0.ID] = AddPenVal * AddViolation1;
                        if (clusterS0.num > UBound)
                        {
                            int AddViolation0 = clusterS0.num - UBound;
                            DropPenalty[clusterS0.ID] = DropInduceVal * AddViolation0;
                        }
                        else
                        {
                            DropPenalty[clusterS0.ID] = 0;
                        }
                    }
                    else if (clusterS0.num <= LBound[clusterS0.ID - 1])
                    {
                        int DropViolation1 = 1 + LBound[clusterS0.ID - 1] - clusterS0.num;
                        DropPenalty[clusterS0.ID] = DropPenVal * DropViolation1;
                        if (clusterS0.num < LBound[clusterS0.ID - 1])
                        {
                            int DropViolation0 = LBound[clusterS0.ID - 1] - clusterS0.num;
                            AddPenalty[clusterS0.ID] = AddInduceVal * DropViolation0;
                        }
                        else
                        {
                            AddPenalty[clusterS0.ID] = 0;
                        }
                    }
                    else
                    {
                        AddPenalty[clusterS0.ID] = 0;
                        DropPenalty[clusterS0.ID] = 0;
                    }
                    TabuEnd[elementT0.ElementID, clusterSA.ID] = Iter + TabuTenure;
                    TabuEnd[elementT0.ElementID, clusterS0.ID] = Iter + SmallTabuTenure;
                    Iter++;
                    #endregion
                    #endregion
                }
                else
                {
                    #region///diversifiction phase
                    MoveLogList.Add("\r\nTS completed. Best solution in Iter No." + (LastBestIter - 1).ToString() + " Start the diversification phase.\r\n************\r\n\r\nStart the Improvement phase.");
                    if (!VarianceConcerned)
                    {
                        VarianceConcerned = true;
                        BestObjVal = BestObjVal * A1 + A2 * BestVariVal;
                        MoveLogList.Add("\r\n\r\n Add Variance. BestObjVal= " + A1.ToString() + "×BestObjVal + " + A2.ToString() + "×Variance = " + BestObjVal.ToString("0.0000") + "\r\n************\r\n\r\n");
                    }
                    int[] tempClusterID = new int[numT];
                    for (int i = 0; i < numT; i++)
                    {
                        Element tempelement = getElementbyID(i + 1);
                        tempClusterID[i] = tempelement.ClusterID;
                        tempelement.MoveTo(this, bestClusterID[i]);
                    }//print the best solution in the last improvement phase
                    AllValuesCalculator();
                    ClusterLogList.Add("\r\n\r\nBest solution in the last improvement phase.\r\nCID\t p\t w\t\t q\t\t w/q\t\t Variance" + '\n');
                    for (int i = 1; i <= numS; i++)
                    {
                        int numbers = getClusterbyID(i).num;
                        int connections = numbers * (numbers - 1) / 2;
                        double weight = Sum_Value[i] / connections;
                        ClusterLogList.Add(i.ToString("") + " \t " + numbers.ToString("") + "\t" + Sum_Value[i].ToString("0.00") + "\t\t" + connections + "\t\t" + weight.ToString("0.0000") + "\t\t" + Vari_Value[i].ToString("0.0000") + '\n');
                    }

                    //rollback
                    for (int i = 0; i < numT; i++)
                    {
                        Element tempelement = getElementbyID(i + 1);
                        tempelement.MoveTo(this, tempClusterID[i]);
                    }
                    //balance the Cluster
                    BalanceCluster();
                    //each Cluster Now has the same number of elements,then change the LBound[] to get a better balanced result
                    for (int i = 0; i < numS; i++)
                    {
                        int newLBound = ((Cluster)getClusterbyID(i + 1)).num / 2;
                        if (LBound[i] < newLBound)
                        {
                            LBound[i] = newLBound;
                        }
                    }
                    InitTSPhase();
                    ClusterLogList.Add("\r\n\r\nClusters after diversifiction phase.\r\nCID\t p\t w\t\t q\t\t w/q\t\t Variance" + '\n');
                    for (int i = 1; i <= numS; i++)
                    {
                        int numbers = getClusterbyID(i).num;
                        int connections = numbers * (numbers - 1) / 2;
                        double weight = Sum_Value[i] / connections;
                        ClusterLogList.Add(i.ToString("") + " \t " + numbers.ToString("") + "\t" + Sum_Value[i].ToString("0.00") + "\t\t" + connections + "\t\t" + weight.ToString("0.0000") + "\t\t" + Vari_Value[i].ToString("0.0000") + '\n');
                    }

                    BestImprove = false;
                    LastImproveIter = Iter;

                    #endregion
                }
            }
            MoveLogList.Add("\r\nThe last TS completed.Best solution in Iter No." + (LastBestIter - 1).ToString());

            #region//roll back to the best solution
            for (int i = 0; i < numT; i++)
            {
                Element tempelement = getElementbyID(i + 1);
                tempelement.MoveTo(this, bestClusterID[i]);
            }
            AllValuesCalculator();
            //recalculate All the values:Val[s] ,ObjVal,AddVal,MoveVal   
            ClusterLogList.Add("\r\n\r\nBest solution in the last improvement phase.\r\nCID\t p\t w\t\t q\t\t w/q\t\t Variance" + '\n');
            for (int i = 1; i <= numS; i++)
            {
                int numbers = getClusterbyID(i).num;
                int connections = numbers * (numbers - 1) / 2;
                double weight = Sum_Value[i] / connections;
                ClusterLogList.Add(i.ToString("") + " \t " + numbers.ToString("") + "\t" + Sum_Value[i].ToString("0.00") + "\t\t" + connections + "\t\t" + weight.ToString("0.0000") + "\t\t" + Vari_Value[i].ToString("0.0000") + '\n');
            }

            #region 正态化逆向
            if (IfNorm)
            {
                Score = (double[,])OriScore.Clone();
                AllValuesCalculator();  
                ClusterLogList.Add("\r\n\r\n******************************************************************************************************"
                    + "\r\nBest solution (Mapped to the initial value).\r\nCID\t p\t w\t\t q\t\t w/q\t\t Variance" + '\n');
                double tempAver = 0;
                double tempVari = 0;
                for (int i = 1; i <= numS; i++)
                {
                    int numbers = getClusterbyID(i).num;
                    int connections = numbers * (numbers - 1) / 2;
                    double weight = Sum_Value[i] / connections;
                    ClusterLogList.Add(i.ToString("") + " \t " + numbers.ToString("") + "\t" + Sum_Value[i].ToString("0.00") + "\t\t" + connections + "\t\t" + weight.ToString("0.0000") + "\t\t" + Vari_Value[i].ToString("0.0000") + '\n');
                    tempAver += weight;
                    tempVari += Vari_Value[i];
                }
                ClusterLogList.Add("Best ObjValue = " + A1.ToString() + "×" + tempAver.ToString("0.0000") +" + "+A2.ToString()+"×"+tempVari.ToString("0.0000")+" = "+(tempAver*A1+tempVari*A2).ToString("0.0000")+ "\n");
            }
            #endregion
            #endregion
        }

        private void BalanceCluster()
        {
            Random ra = new Random();
            int basenum = numT / numS;
            int[] temp = new int[numS];
            for (int i = 0; i < numS; i++)
            {
                Cluster tempCluster = getClusterbyID(i + 1);
                temp[i] = tempCluster.num - basenum;
            }
            //move extra element to dummy Cluster
            for (int i = 0; i < numS; i++)
            {
                Cluster curCluster = getClusterbyID(i + 1);
                while (temp[i] > 0)
                {
                    int randomOrder = ra.Next(1, curCluster.num + 1);
                    Element curElement;
                    curElement = curCluster.First;
                    for (int j = 0; j < randomOrder - 1; j++)
                    {
                        curElement = curElement.After;
                    }
                    curElement.MoveTo(this, 0);
                    temp[i]--;
                }
            }
            //move elements from dummyCluster to the cluster in short 
            Cluster dumCluster = getClusterbyID(0);
            Element dumElement = null;
            for (int i = 0; i < numS; i++)
            {
                //Cluster curCluster = getClusterbyID(i + 1);
                while (temp[i] < 0)
                {
                    dumElement = dumCluster.First;
                    dumElement.MoveTo(this, i + 1);
                    temp[i]++;
                }
            }
            dumElement = dumCluster.First;
            while (dumElement != null)
            {
                dumElement.MoveTo(this, ra.Next(1, numS + 1));
                dumElement = dumCluster.First;
            }
        }


        private string BugCheck3()
        {
            double tempOBJVal = 0;
            string slog = "";
            for (int s = 1; s <= numS; s++)
            {
                Cluster curcluster = getClusterbyID(s);
                if (curcluster.NumberLinks != (curcluster.num * (curcluster.num - 1) / 2))
                {
                    slog = slog + "\tERROR:NumLinks in Cluster(" + s.ToString() + ") is " + curcluster.NumberLinks.ToString() + " ,which should be " + (curcluster.num * (curcluster.num - 1) / 2).ToString() + "\r\n";
                }
                double temp = Calculator(curcluster);
                if (Sum_Value[s] != temp)
                {
                    slog = slog + "\tERROR:Sum_Value(" + s.ToString() + ") is " + Sum_Value[s].ToString() + " ,which should be " + temp.ToString() + "\r\n";
                }
                //Cluster's numberlinks&value have been comfirmed above
                tempOBJVal += Sum_Value[s] / curcluster.NumberLinks;
                for (int t = 1; t <= numT; t++)
                {
                    Element curelement = getElementbyID(t);
                    temp = AddValueCalculator(curelement, curcluster);
                    if (AddVal[t, s] != temp)
                    {
                        slog = slog + "\tERROR:AddValue(" + t.ToString() + " , " + s.ToString() + ") is " + AddVal[t, s].ToString() + " ,which should be " + temp + "\r\n";
                    }
                    temp = MoveValueCalculator(curelement, curcluster);
                    if (MoveVal[t, s] != temp && MoveVal[t, s] != DropPenVal)
                    {
                        slog = slog + "\tERROR:MoveValue(" + t.ToString() + " , " + s.ToString() + ") is " + MoveVal[t, s].ToString() + " ,which should be " + temp + "\r\n";
                    }
                }
            }
            if (!MathHelper.IfnearlyEqual(tempOBJVal, ObjVal))
            {
                slog = slog + "\tERROR:ObjVal is " + ObjVal.ToString() + " ,which should be " + tempOBJVal.ToString() + "\r\n";
            }
            return slog;
        }


        public Cluster getClusterbyID(int id)
        {
            return (Cluster)clusterList[id];
        }
        public Element getElementbyID(int id)
        {
            return (Element)elementList[id - 1];
        }


        public double CalculateScore(Element t1, Element t2)//坐标距离计算  当前版本不需要计算，直接从Score中取值
        {
            //double score = 0;
            //if (t1 == t2)
            //{
            //    return 0;
            //}
            //if (t1 == null || t2 == null)
            //{
            //    return double.MaxValue;
            //}
            //for (int i = 0; i < t1.AttributesX.Count; i++)
            //{
            //    double tempscore = Convert.ToDouble(t1.AttributesX[i]) - Convert.ToDouble(t2.AttributesX[i]);
            //    if (tempscore < 0)
            //    {
            //        tempscore = 0 - tempscore;
            //    }
            //    score +=  tempscore*tempscore;
            //}
            //int X = t1.ElementID;
            //int Y = t2.ElementID;
            //score = Score[X, Y];
            //if (score > 0)
            //{
                //return score;
            //}
            //else
            //{
            //    return -score;
            //}
            return Score[t1.ElementID, t2.ElementID];
        }

        public double Calculator(Cluster c)//return Value of Cluster s
        {
            double temp = 0;
            Element em1 = c.First;
            Element em2;
            while (em1 != null)
            {
                em2 = em1.After;
                while (em2 != null)
                {
                    temp += CalculateScore(em1, em2);
                    em2 = em2.After;
                }
                em1 = em1.After;
            }

            return temp;

        }

        private void AllValuesCalculator()
        {
            ObjVal = 0;
            for (int s = 1; s <= numS; s++)
            {
                Cluster curcluster = getClusterbyID(s);
                curcluster.NumberLinks = (curcluster.num * (curcluster.num - 1) / 2);
                double temp = Calculator(curcluster);
                Sum_Value[s] = temp;
                Vari_Value[s] = VariValuesCalculator(curcluster);

                //Cluster's numberlinks&value have been comfirmed above
                if (VarianceConcerned)
                {

                    ObjVal += A1 * (Sum_Value[s] / curcluster.NumberLinks) + A2 * Vari_Value[s];
                }
                else
                {
                    ObjVal += Sum_Value[s] / curcluster.NumberLinks;
                }
                for (int t = 1; t <= numT; t++)
                {
                    Element curelement = getElementbyID(t);
                    temp = AddValueCalculator(curelement, curcluster);
                    AddVal[t, s] = temp;
                    temp = MoveValueCalculator(curelement, curcluster);
                    MoveVal[t, s] = temp;
                    temp = Move_Vari_ValueCalculator(curelement, curcluster);
                    Move_Vari_Val[t, s] = temp;
                }
            }
        }

        private double FullvaluesCalculator(Cluster c)
        {
            double values = Calculator(c);
            double fullvalues = values / c.NumberLinks;
            return fullvalues;
        }

        private double VariValuesCalculator(Cluster c)//should confirm SumValue updated!!
        {
            double fullvalue = Sum_Value[c.ID]/c.NumberLinks;
            double varivalue = 0;
            Element ele1, ele2;
            ele1 = c.First;
            while (ele1 != null)
            {
                ele2 = ele1.After;
                while (ele2 != null)
                {
                    double temval = CalculateScore(ele1, ele2);
                    varivalue += (temval - fullvalue) * (temval - fullvalue);
                    ele2 = ele2.After;
                }
                ele1 = ele1.After;
            }
            return varivalue / c.NumberLinks;
        }

        private double ObjValueCalculator()
        {
            double objvalues = 0;
            for (int i = 1; i < clusterList.Count; i++)
            {
                double tempfullvalues = FullvaluesCalculator(getClusterbyID(i));
                objvalues = objvalues + tempfullvalues;
            }
            return objvalues;
        }

        private double AddValueCalculator(Element t0, Cluster s0)
        {
            Element t = s0.First;
            double tempval = 0;
            while (t != null)
            {
                tempval += CalculateScore(t0, t);
                t = t.After;
            }
            return tempval;
        }

        private double MoveValueCalculator(Element t0, Cluster sA)
        {
            double valueSA = Calculator(sA);
            double _valueSA = valueSA + AddValueCalculator(t0, sA);
            int numlinkSA = sA.num * (sA.num - 1) / 2;
            int _numlinkSA = numlinkSA + sA.num;
            double DelFullA = _valueSA / _numlinkSA - valueSA / numlinkSA;

            Cluster s0 = getClusterbyID(t0.ClusterID);
            double valueS0 = Calculator(s0);
            double _valueS0 = valueS0 - AddValueCalculator(t0, s0);
            int numlinkS0 = s0.num * (s0.num - 1) / 2;
            int _numlinkS0 = numlinkS0 - s0.num + 1;
            double DelFullD = _valueS0 / _numlinkS0 - valueS0 / numlinkS0;

            return DelFullA + DelFullD;
        }


        private double Move_Vari_ValueCalculator(Element t0, Cluster sA)//should confirm AddVal[t0,s0]&AddVal[t0,sA] updated!!
        {
            Cluster s0 = getClusterbyID(t0.ClusterID);

            if (s0.ID == sA.ID)
            {
                return 0;
            }
            double NewFullValA = (Sum_Value[sA.ID] + AddVal[t0.ElementID, sA.ID]) / (sA.NumberLinks + sA.num);
            double NewFullVal0 = (Sum_Value[s0.ID] - AddVal[t0.ElementID, s0.ID]) / (s0.NumberLinks - s0.num + 1);

            double new_Vari_s0 = 0;
            Element ele1, ele2;
            double temval;
            #region MoveVariValue Changed in S0
            ele1 = s0.First;
            while (ele1 != null)
            {
                if (ele1 != t0)
                {
                    ele2 = ele1.After;
                    while (ele2 != null)
                    {
                        if (ele2 != t0)
                        {
                            temval = CalculateScore(ele1, ele2);
                            new_Vari_s0 += (temval - NewFullVal0) * (temval - NewFullVal0);
                        }
                        ele2 = ele2.After;
                    }
                }
                ele1 = ele1.After;

            }
            #endregion

            #region MoveVariValue Changed in SA
            double new_Vari_SA = 0;
            ele1 = sA.First;
            while (ele1 != null)
            {
                ele2 = ele1.After;
                while (ele2 != null)
                {
                    temval = CalculateScore(ele1, ele2);
                    new_Vari_SA += (temval - NewFullValA) * (temval - NewFullValA);
                    ele2 = ele2.After;
                }
                temval = CalculateScore(ele1, t0);
                new_Vari_SA += (temval - NewFullValA) * (temval - NewFullValA);
                ele1 = ele1.After;
            }

            #endregion
            new_Vari_s0 = new_Vari_s0 / (s0.NumberLinks - s0.num + 1);
            new_Vari_SA = new_Vari_SA / (sA.NumberLinks + sA.num);


            return new_Vari_s0 + new_Vari_SA - Vari_Value[s0.ID] - Vari_Value[sA.ID];
        }

        public ArrayList clusterList;
        public ArrayList elementList;
        //private double StartWeight;//初始权重(为凸显某个因素的权)设为0.7

        private ArrayList[] CandidateList;
        private ArrayList MoveLogList;
        public ArrayList GetMoveLogList() { return MoveLogList; }
        private ArrayList ClusterLogList;
        public ArrayList GetClusterLogList() { return ClusterLogList; }
        private ArrayList ObjValList, Obj_Vari_ValList;
        public ArrayList GetObjValList() { return ObjValList; }
        public ArrayList GetObj_Vari_ValList() { return Obj_Vari_ValList; }

        //private double Weight;
        private double ObjVal;
        private double Obj_Vari_Val;

        //attributes for TS(penalties/positive) 
        private double DropPenVal, AddPenVal, AddInduceVal, DropInduceVal;
        private int NumUBViolate, NumLBViolate, NumAllViolate;
        private bool BestImprove, VarianceConcerned;

        private double[,] TabuEnd;
        private double BestObjVal, BestVariVal;
        private double SmallTabuTenure, TabuTenure;
        private int LastImproveIter, Iter, FirstBestIter, LastBestIter, NumNewBest;

        private int numS;
        private int numT;
        private int UBound;
        private int[] LBound;
        private double A1, A2;//weight of Average&Variance

        private double[,] MoveVal;
        private double[,] Move_Vari_Val;
        private double[,] AddVal;

        public double[] AddPenalty;//trialVer
        public double[] DropPenalty;//trialVer
        public double[] Sum_Value;//trialVer
        public double[] Vari_Value;//trialVer


        private double[,] Score;//trialVer
        private double[,] OriScore;//trialVer
        bool IfNorm;
        public void ImportScore(double[,] score, bool ifnormalized, int ElementNum)
        {

            Score = (double[,])score.Clone();
            OriScore = (double[,])score.Clone();
            if (ifnormalized)
            {
                IfNorm = true;
                MathHelper.MaxMinNorm(Score, ElementNum + 1, ElementNum + 1);
            }
        }
        public void ImportScore(int[,] score)
        {
            Score = (double[,])score.Clone();
        }
    }

}
